library(ggplot2)
library(dplyr)
library(tidyr)
library(ggrepel)
library(multiplot)
#染色质状态
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201222")
rt=read.csv("HMM.enrich.117K.rm.eQTL.support.psyASH.meta.csv",header=T)
rt$group=ifelse(rt$P.value<0.05,"sig","nosig")
rt1=rt[rt$NUM.state<10,]
rt2=rt[rt$NUM.state>=10,]
rt1=arrange(rt1,group,P.value)
rt2=arrange(rt2,group,P.value)
rt=rbind(rt2,rt1)
rt$num=1:15
rt$logOR=log2(rt$OR)
rt$logupper=log2(rt$upper)
rt$loglower=log2(rt$lower)


p1=ggplot(data=rt,aes(y=logOR,x=num))+geom_errorbar(aes(ymin=loglower,ymax=logupper,group=group,color=group),width=0.5)+geom_point(aes(size=0.5,color=group))+
    scale_color_manual(values = alpha(c("#3C5488","#E64B35"),0.9))+coord_flip()+geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_x_continuous(breaks = rt$num,labels = rt$DESCRIPTION)+
    theme_bw()+ylab("log2(OR)")+
    theme(panel.grid.minor = element_blank(),panel.grid.major = element_blank(),axis.text.x = element_text(size = 15,color="black"),axis.text.y = element_text(size = 15,color="black"))

#H3k的组蛋白甲基化修饰出现的可能性与羟甲基化方向，属于调控类型
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201227.H3k.analysis/807.analysis/")
file=read.csv("807.ASH.analysis.csv",header =T)
file$group=ifelse(file$P.value<0.05,"significant","nonsignificant")
file=file[order(file$P.value,decreasing = T),]
file1=data.frame(file[6:7,1:3])
file1$num=1:2
file1=gather(file1,key = direction.group,value = counts,same.mean,opposite.mean)

p2=ggplot(file1,aes(x=num,y=counts,,group=direction.group))+geom_bar(stat ="identity",width = 0.5,position = "dodge",aes(fill=direction.group))+
  scale_fill_manual(values = alpha(c("#3C5488","#E64B35"),1))+scale_x_continuous(breaks = file1$num,labels = file1$H3k.group)+	
  geom_text(aes(label = counts),position=position_dodge(width = 1),size = 5,vjust = -0.25)+coord_flip()+theme_classic(base_size = 15)+labs(x="",y="counts",title = "")+guides(fill=F)#这里去掉图例是为了更好地画图，后期会在Ai里复制一个图例过去

###TF analysis 采用807预测的TF，然后分析其与羟甲基化的拟合关系，用于判断羟甲基化状态是否会影响TF的亲和力。
library(ggplot2)
library(ggrepel)
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20210103.motif.analysis")
#1.全部TF做散点图
rt=read.csv("cor.mean.vaf.vs.score.allele.diff.csv",header = T)
rt$group="no.prefer"
sel1=which(rt$cor.p.value<0.05&rt$cor.estimate>0.7)
sel2=which(rt$cor.p.value<0.05&rt$cor.estimate< -0.7)

rt[sel1,]$group="prefer.alt.allele"
rt[sel2,]$group="prefer.ref.allele"
rt$logpvalue=-log10(rt$cor.p.value)
rt$absR=abs(rt$cor.estimate)

p3=ggplot(rt,aes(y=absR,x=logpvalue,color=group))+geom_vline(xintercept = -log10(0.05),color="red",linetype="dashed")+
    geom_point(size=3)+scale_color_manual(values = alpha(c("#666666","#3C5488","#E64B35"),0.9))+
    labs(y="r",x="-log10 (Pvalue)",title="")+
    theme_classic(base_size = 15)+	theme(legend.position = "none")+	#这里去掉图例是为了更好地画图，后期会在Ai里复制一个图例过去
    geom_text_repel(data=subset(rt,rt$cor.p.value<0.05),aes(label=TF),size=3, fontface="bold",force = T, segment.color = "black", show.legend = FALSE)

  cor_eqn=function(data){
    r=cor.test(data$mean.vaf,data$score.allele.diff,method = "pearson")
    eq <- substitute(italic(P.value) == p ~~ ","~~italic(r)~"="~rr,
                     list(p=as.numeric(format(r$p.value, digits = 3)),
                          rr=as.numeric(format(r$estimate, digits = 3))
                     ))
    as.expression(eq)
}					#粘贴在图中的公式

result=read.table("807.psy.ASH.motif.VAF.txt",head=T,sep="\t")
result$id2=paste(result$rsid,result$geneSymbol,sep="_")
result=result[!duplicated(result$id2),]
result$score.allele.diff=abs(result$scoreAlt-result$scoreRef)
tmp=result[result$geneSymbol=="NRF1",]

r=cor.test(tmp$mean.vaf,tmp$score.allele.diff,method = "pearson")
p=as.numeric(format(r$p.value, digits = 3))
rr=as.numeric(format(r$estimate, digits = 3))

p4=ggplot(tmp,aes(y=score.allele.diff,x=mean.vaf))+geom_point(size=3)+
  geom_smooth(method = "lm",color="black")+labs(y="binding likehood differences",x="VAF",title="")+
  theme_classic(base_size = 15)+annotate(x=0.35,y=9.9,label=cor_eqn(tmp),geom ="text" )

dt=data.frame(table(rt$group))
dt[dt$Var1=="no.prefer",]$Freq=183+dt[dt$Var1=="no.prefer",]$Freq

p5=ggplot(dt, aes(x = "", y = Freq, fill = Var1)) +	#这是普通的饼图
  geom_bar(stat = "identity", width = 1)+
  coord_polar(theta = "y") + 
  labs(x = "", y = "", title = "") + 
  theme(axis.ticks = element_blank()) + 
  theme(legend.title = element_blank(), legend.position = "top") + 
  scale_fill_discrete(breaks = dt$Var1, labels = dt$Freq) + scale_fill_manual(values = alpha(c("#666666","#3C5488","#E64B35"),0.9))+theme_bw()+
  theme(
    panel.grid = element_blank(),
    panel.border= element_blank(),
    axis.text.y = element_blank(),
    axis.text.x = element_blank(), 
    axis.ticks = element_blank(),
    axis.title = element_blank()
  )+guides(fill=F) 
  
  result=dt
result$ratio=round(result$Freq/sum(result$Freq),4)
result=result[order(result$Freq),]
result$lognum=log(result$ratio*1000)	#玫瑰图
p5=ggplot(result,aes(x=reorder(Var1,Freq), y=lognum, fill=Var1))+ geom_bar(stat="identity", color="black")+coord_polar()+theme_bw()+
  ylim(-5,12)+theme(
    panel.grid = element_blank(),
    panel.border= element_blank(),
    axis.text.y = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks = element_blank(),
    axis.title = element_blank()
  )+
  geom_text(aes(label = paste0(Var1,"\n",Freq)),data = result,vjust = "left", hjust = "outward", color = "dark blue", fontface="bold")+
  scale_fill_manual(values = alpha(c("#666666","#3C5488","#E64B35"),0.9))+guides(fill=F)


  
    layout <- matrix(c(1,1,2,3,1,1,2,3,1,1,5,5,4,6,5,5,4,6,5,5), nrow = 4)
multiplot(plotlist=list(p1,p5,p4,p2,p3),layout=layout)
